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1.  INTRODUCTION 


In  the  development  of  the  computational  tools  to  study  radio  wave  prop¬ 
agation  in  the  ionosphere,  it  is  useful  to  have  exact  ray  calculations  through 
realistic  electron  density  profiles.  The  parabolic  model  ionospheric  layer 

was  introduced  by  Fosterling  and  Lasser^^  in  1931,  and  further  developed 

(2)  (3) 

and  refined  by  Appleton  and  Beynon,  '  and  Rawer.  '  This  model  leads  to 
exact  solutions  for  the  vertical  ionogram  and  for  horizontally  stratified 
layers,  but  requires  an  iterative  technique  for  spherical  geometry. 

(4) 

A.  H.  DeVoogt  described  the  quasi-parabolic  layer  as  an  alternative 
computational  tool,  which  was  developed  in  detail  by  T.  A.  Croft  and 
H.  Hougasian!^  This  model  layer  leads  to  representation  of  the  solut¬ 
ions  to  the  ray  equations  in  terms  of  the  elementary  functions. 

In  addition,  D.  B.  Muldrew^  introduced  a  quasi  linear  election  density 

(7) 

profile,  which  was  developed  in  detail  by  D.  E.  Westover. 

/  o\ 

P.  A.  Bradley  and  J.  R.  Dudeney  '  developed  a  parabolic/linear 

complete  election  density  profile  as  a  candidate  for  a  standard  profile. 

(9) 

J.  D.  Milsom  modified  the  Bradley-Dudeney  profile  to  be  represented 
by  Croft  and  Hougasian  type  quasi  parabolic  segments  connected  by  a 
Westover  type  quasi  linear  layer. 

All  these  models  are  strictly  valid  for  ionospheres  which  display  spher¬ 
ical  symmetry.  A  major  problem  is  to  treat  in  a  systematic  way  the  pres¬ 
ence  of  horizontal  gradients  in  the  electron  density.  K.  Folkestad^*^ 
described  a  method  of  correcting  the  results  of  ray  calculations  through 
the  assumption  of  tilts  in  the  ionosphere.  This  technique  assumes  a  spher¬ 
ically  symmetric  ionosphere,  except  that  the  center  of  curvature  is  assumed 


not  necessarily  to  be  at  the  center  of  the  earth. 

The  problem  of  determining  the  center  of  curvature  of  the  ionosphere 
was  addressed  by  R.  I.  Beckwith^  ^  who  used  the  spatial  variation  of 
virtual  height  with  parabolic  layer  theory  calculations. 

D.  Odom  (unpublished)  used  the  election  density  gradient  in  the 
neighborhood  of  the  reflection  point  to  construct  approximate  corrections 
to  the  ray  trajectories  for  the  effects  of  horizontal  gradients. 

No  simple  technique  appears  to  be  forthcoming  to  handle  the  effects 

of  the  geomagnetic  field  or  neutral  collisions;  however,  ionospheric 

absorption  may  be  calculated  for  field  strength  calculations  by  using  a 

( 1 2) 

semi-empirical  model  such  as  the  Bradley-George  '  absorption  model. 

The  next  section  is  devoted  to  a  development  of  the  ray  equations  from 
the  Hamiltonian  equations  described  by  J.  Haselgrove-13^  the  approach  is, 
firstly,  ignoring  the  magnetic  field  and  collisions,  and,  secondly,  from  an 
assumption  of  nearly  spherically  symmetric  ionosphere.  Equations  are 
developed  which  allow  an  approximate  correction  for  the  effects  of  hori¬ 
zontal  gradients,  which  may  be  applied  after  calculation  of  the  major  ray 
parameters  through  a  given  layer. 

Section  3  describes  the  quasi-parabolic  ray  calculation  equations, 
and  the  options  available  for  segmenting  profiles  with  quasi-parabolic, 
inverted  quasi  parabolic,  and  quasi-linear  layers. 

Section  4  describes  the  Bradley-Dudeney  profile  and  suggested  modi¬ 
fications  to  quasi  parabolic/quasi  linear  models,  and  Section  5  describes 
the  piecewise  smooth  quasi  parabolic  profile  with  various  construction 
options. 


In  Section  6  we  discuss  models  for  simulation  of  ionospheric  para¬ 
meters,  updating  techniques  based  on  real  time  ionograms,  and  profile 
modifications  to  account  for  horizontal  gradients.  Section  7  is  devoted 
to  applications  of  the  model  to  transmission  circuits,  while  Section  8 
describes  applications  to  vertical  and  oblique  ionogram  simulation. 

Section  9  is  devoted  to  program  descriptions  and  listings  and  may 
serve  as  a  user  guide  to  the  program  modules  implementing  the  ray 
calculations. 


2.  DEVELOPMENT  OF  THE  RAY  EQUATIONS 


When  the  effects  of  the  earth's  magnetic  field  and  electron  collisions 

/  j  2) 

may  be  neglected,  the  equations  governing  electromagnetic  wave  tra¬ 

jectories  are  given  by. 


5  =  V6/r 
CP  =  V  / r  sine 

V  =  n  —  +  VQ  e  +  V„sin0m 

r  ^  ar  e  cp  ^ 

rVe  =  M  M  -  V  +  rVcpcos9£e 

r  sineV^  =  p  -  Vising  r  -  rV^cos^e 


(2.  1) 
(2.2) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 


where  r,  9,  cp  are  the  spherical  polar  coordinates  of  a  point  on  the  ray  path, 
Vr,  Vg,  are  the  components  of  the  ray  normal  direction  in  the  r,  9, 
and  Cp  directions  so  normalized  that 


Vr  +  V9  +  M  <2'7» 

The  phase  refractive  index  p,  in  the  no  field,  no  collision  case  is  given  by 


2 

M 


1  -  (f 


(2.8) 


where  f  is  the  local  plasr  a  frequency  and  f  is  the  radiation  frequency;  in 
'  P 

in  this  case,  the  group  refractive  index  n'  =  ^j-(fM)  is  given  by 


M1  =  1/M 


(2.9) 


The  independent  variable  t  is  the  group  path  defined  by  (no  field,  no  collision) 


t  =  „  p'dg 


(2.  10) 
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Equation  (2.6)  may  be  immediately  rewritten  as 


sinO 


V 


1  Sii2 

2  Sep 


and  equations  (2. 


where 

and 


5)  and  (2,  6)  may  be  combined  as 


-4-(rV  )2 

dt  '  FT 


r2  e 


5m 


dep  5m 

"3T  W 


TV 

cp 


VH  -  Ag  +  Vl  =  H  cosS 
Pis  the  local  elevation  angle. 


(2.  11) 


(2.  12) 


In  the  absence  of  horizontal  gradients  in  the  refractive  index, 

Equation  (2.  11)  expresses  the  requirement  of  great  circle  propagation, 

and  Equation  (2.  12)  is  the  well  known  Bouger's  rule,  pr  cos  3  is  a 

constant  of  the  motion.  Furthermore,  the  ray  trajectory  equations 

can  be  reduced  to  quadratures  in  the  radial  coordinate, 

r 

0(r)  =  (m  r  cos  g  )  0  dr' 

°  0  0  ra  r'  Z(r ’) 


(2.  13) 


P(r) 


r"  r'  dr' 


Jr  Z<r’) 


(2.  14) 


P(r) 
Z  (r) 


2 

p  r1  M  (r1)  dr1 

°ro  Z<r’> 

M (r)  r  sin  g  (r) 


(2. 15) 

2]  1/2 

(M  r  cof  g  )  (2.  16) 

0  0  o 


where  0(r),  P'(r),  and  P(r)  are  respectively  the  ground  range  (in  radians), 
the  group  path  length,  and  the  phase  path  length  accrued  along  the  ray 
trajectory  between  entering  a  layer  at  r  and  the  point  on  the  trajectory 

O 

at  height  r.  Integration  extends  to  ray  apogee  or  perigee  (Z  =  0)  or  to 


penetration  through  the  layer,  whichever  is  appropriate. 

If  the  horizontal  gradients  in  the  index  of  refraction  are  small  enough, 
the  Equations  (2.  13)  through  (2.  15),  suitably  implemented,  remain  good 
approximations  to  the  ray  path  parameters.  The  effects  of  horizontal 
components  of  gradients,  then,  may  be  estimated  and  applied  as  cor¬ 
rections  to  the  local  elevation  angle  and  deviation  to  great  circle  prop¬ 
agation.  From  Equation  (2.8),  partial  derivatives  of  the  phase  refractive 
index  may  be  written  as 

=  2(p2-  1)-J-  (log  f  )  (2.17) 

ax  '  &x  p 

Consider  first  the  case  where  the  component  transverse  to  the  plane 
of  ray  propagation  vanishes,  and  let  9  be  the  ground  range  and  cp  the  rela¬ 
tive  bearing  from  that  plane.  Great  Circle  propagation  still  pertains, 
thus,  VCp(t)  =  0,  cp(t)  =  0,  vH(t)  =  r®,  and  from  Equation  (2.  12), 

■^-(pr  cosP)  =  -^-(l°g  ^p)(*J  "  1)  (2.  18) 

If  the  longitudinal  gradient  of  log  f  can  be  considered  sufficiently  constant, 
then  the  above  may  be  integrated  to  give  the  change  in  the  Bouger  invariant 
due  to  this  gradient  accrued  through  the  layer. 

A(pr  cosf3)  =  (P  -  P')  (log  f  )  (2.  19) 

This  correction,  then,  should  be  ft  1  back  into  the  integrals.  Equations  (2.  13) 
through  (2.  15),  for  recalculation.  This  process  may  be  reiterated,  if  desired, 
until  stability  of  the  computation  is  achieved. 

In  the  presence  of  transverse  gradients  Equation  (2.  11)  can  be  similarly 
integrated  to  provide  the  accrual  to  the  invariant  r  sin  thus, 

A(r  sin9  V  )  -  (P  -  P«)  (log  f  )  (2.20) 

w  o  w  p 


9 


or 


r  sin  6  V  =  r  sin©  V  +  (P  -  P’)-r^<log  f  ) 

CD  o  o  tPo  ^  V  P 


from  which, 


r  sin 6  Vcp 
]Jr  cosp 


] /sin  6 


where 


cp(  e  )  =  CP  4  A  Kjn(  0  -  9) 
c  o 


r  s  i  n  6  V  cp 


(2.21) 


(2.22) 


M  r  cosg  AVG  sinS  sinS 


(2.23) 


Finally  the  contribution  from  the  transverse  gradient  can  be  included 


in  Equation  (2.  19),  thus 


MMr  cosg)  =  (P  -  P’){ye-l°g  fp  4  A  lVloe  fP} 


(2.24) 


Clearly,  any  implementation  of  these  last  few  equations  in  numerical 
calculations  must  be  careful  to  avoid  the  kinematic  singularity  at  0=  0. 
The  best  technique  would  be  to  arrange  the  calculation  so  that  the  reflect¬ 
ion  point  occurs  approximately  tt/2  from  the  pole  of  the  coordinate 
system. 


10 
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3.  INTEGRATION  OF  THE  RAY  EQUATIONS 
THROUGH  QUASI-PARABOLIC  LAYERS 


The  parabolic  model  of  election  density  has  been  for  many  years  very 

popular  as  a  technique  to  simulate  and  study  the  phenomenom  of  ionospheric 

2 

radio  propagation.  This  model  expresses  the  electron  density  kf^  by  the 
equation  of  a  parabola,  to  wit, 


fn  (r)  =  fc 
P  c 

=  0 


r  -y  <r  r 

m  **  m 

elsewhere 


+  y, 


m 


(3.  1) 


where  f^(r)  is  the  plasma  frequency  having  a  maximum  value  of  f^  at  the 
layer  semi -thickness,  and  r  is  the  radial  distance  from  the  center  of  the 
earth. 


For  vertical  incidence  (yir  cos  3=0)  this  model  allows  solutions  for 
0,  P',  and  P  in  terms  of  elementary  functions,  for  in  this  case 
Z(r)  =  r  M (r) 

and  the  integrals  for  P'  and  P  (©  =  0  trivially)  take  the  form 


P'(r)  = 


(i-T* - > 

JVa  r  +  br  +  c 

f  I z - 1 

P  (r)  =  I  y  a  r  +br  +  c  dr 


(3.2) 


However,  for  oblique  propagation  the  expression  under  the  radical  is  a 
quartic,  and  S,  P',  and  P  are  represented  by  elliptic  integrals  which  do 
not  allow  for  representation  in  terms  of  the  elementary  functions.  On  the 
other  hand,  the  virtual  height  method  has  been  developed  to  approximate 
the  solution  for  P'  on  the  oblique  path  from  the  vertical  solution  for  the 
equivalent  vertical  frequency. 


As  pointed  out  by  several  investigators  (References  4,  5,  6,  7,  9)  a 
minor  modification  to  the  parabolic  model  permits  the  representation  of 


the  ray-path  parameters  in  terms  of  elementary  functions.  This  modifi¬ 
cation  is  known  as  the  Quasi- Parabolic  (QP)  Model  layer  and  given  by 


=  0 


rb<  r<rb 


elsewhere 


(3.3) 


where  r,  =  r  -  y 

b  m  m 

and  r '  =  r  +  y 
b  m  m 


r  /it  - 
m  /  \  m  m 


)] 


For  typical  values  of  r  and  y  (6800  and  100  respectively)  the  additional 


m 


m 


term  r^/r  would  cause  a  difference  of  a  few  percent  at  most  from  the  para¬ 
bolic  layer,  and  primarily  in  the  topside  of  the  layer. 

The  advantage  of  the  QP  model  is  clearly  seen  if  we  note  that 


Z^(r)  =  (M  r)^  -  (p  r  cos  (3  )^  =  rZ  (A*  ^  -  2  Bx  +  C) 
o  o  o  m 


(3.4) 


where 


r/r 


m 


A  =  1  •  XC  +  XG 
B  =  X_ 


C  =  X, 


XC=  <fc/f) 

XC=  XC  (rm^m  *  1)2 

T  =  (p  r  cos  P  /r  )‘ 
o  o  o  m 


The  equations  for  ground  range,  group  path  and  phase  path,  then,  have 
solutions  in  terms  of  elementary  functions  as  follows 


®(r)  =  r  F. 


P 
P  ( 


,,r»'  [F.  tBF3]  IA 

r)-  'JFrBF,(  ic  +  r,F2j 


(3.  b) 


)  2 


where 


Fj  -  Y(x)  -  Y(xq) 


F,  =  fx  dx'/x'  Y(x«) 
L  J  x 


Y (x)  =  (Ax2  -  2Bx  +  C)1  ,Z 


(3.6) 


These  may  be  integrated  to  yield 

jx<Bxo  -  C 


J-LOG  - — nr— 


c  -  yY) 


yY)  j 


=  (x  Y  -  xY  )/(Bxx  ) 
'  o  o  o 


=  (x  -  x  )/axx„ 
O  o 


,  |  Y  Y (Bx  -C)  -  Y  Y  (Bx-C) 

=  J-ATAN  < - °  ° 

V  j(Bxo 


C)  (Bx  -  C)  -  CYq  Yj 


—  LOG  < 
a  \ 


B  -  Ax  -  aY 

B  -  Ax  -  aY  i 
o  o 


(Yo  -  Y)/B 


(x  -  xq)/y 


•ATAN 


aY0  (B  -  Ax)  -  aY  (B-AxJ 
(B  -  Ax)  (B  -  Axq)  -  A  Y  Y 


S 


where 

x 

=  r/r  ,  x  =  r  /r 

m  o  o  m 

and 

2 

a 

• 

U_ 

ii 

•V 

ii 

C>0 


c=o 

Bj/0 


C  =  0 
B=0 


C<0  (3.7) 


A>  0 


A=  0 
0 


A=  0 
B=  0 


A<0 

(3.8) 


The  ray  is  assumed  to  enter  the  ionospheric  layer  at  r  ,  and  the  upper 
limit  r  is  the  other  boundary  of  the  layer,  unless  reflection  occurs.  This 
condition  occurs  when  Z  vanishes  (Z  =  p  r  sin(3) 

or  Ax2  -  2B  x^  +  C  -  0  (3.9) 


13 


If  both  solutions  to  this  quadratic  equation  are  physically  valid,  then  one 
solution  belongs  to  the  ray  entering  from  below  and  reflecting  downwards, 
while  the  other  belongs  to  the  ray  entering  from  above  and  reflecting  up¬ 
wards.  After  reflection  the  ray  would  be  carried  back  to  the  appropriate 
boundary  by  doubling  the  accumulated  group  path,  phase  path  and  range. 

When  a  QP  layer  is  employed  as  a  model  for  the  E  or  F  ionospheric 

layers  (y  ^  r  ,F  ^  F  )  A  is  always  positive  while  C  remains  posit- 
7  7  m  <  <  m  p^c  7  r  r 

ive  up  to  the  UHF  band,  thus  the  LOG  forms  prevail  for  functions  F^  and  F^. 
On  the  other  hand  one  may  desire  to  turn  the  parabola  over  for  use  as 
bottomside  or  topside  tails,  or  as  a  fitting  layer  between  two  layers  with 
local  maxima.  This  inverted  QP  layer  is  given  by 

X(r)  =  (F  /f)2  =  X-  +  X  (t/t  -l)2  (3.  1 

p  L,  A  m 

with  X^  2  0  X^  >0. 

In  this  case. 


A 

B  =  -  X 


1  -  X.  -  X 
L  A 


c  =  -XA  -r 


(3.  1 


Clearly  C  is  always  negative  while  A  is  negative  except  at  high  frequencies 
thus  the  ATAN  forms  prevail  for  functions  F^  an<3  F^. 

Another  model  layer  with  this  same  class  of  solutions  is  the  quasi- 
linear  layer  given  by 

X  (r)  -  X_  4  2  XT  (1  -  ,/r) 


(3.12 


in  which  case 


A  =  (1  -  Xm  ;  2  XL) 

B  =  ±  2  XL 

c  =  -  r 

Clearly  C  is  negative,  thus  the  ATAN  form  applies  for  function  F£  while 
A  may  be  positive  or^  negative.  Note  that  the  quasi-linear  layer  defined 

(7) 

by  Equation  (3.  12)  is  not  the  same  as  that  treated  by  Westover  and  by 
Milsom/^ 

These  three  functional  forms,  the  quasi-parabolic,  the  inverted  quasi¬ 
parabolic,  and  the  quasi-linear,  form  a  stable  of  candidates  for  segments 
of  ionospheric  layers  which  allow  the  construction  of  election  density  pro¬ 
files  with  substantial  variation;  one  may  construct,  segment  by  segment, 
peaks  and  valleys  and  tails.  For  any  such  profile  the  contributions  to 
ground  range,  group  path,  and  phase  path  may  be  calculated  in  each 
segment,  from  essentially  the  same  class  of  solutions,  and  then  summed 
for  the  entire  ray  trajectory. 
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4.  THE  MODIFIED  BRADLEY  DUDENEY  PROFILE 

P.  A.  Bradley  and  J.  R.  Dudeney^  described  an  idealized  model  for 
election  density  profiles  generated  from  four  ionospheric  parameters, 
f  E,  f  F~,  h  F_  and  y  F_  (see  Table  1  for  definitions  of  ionospheric, 
parameter  notation).  This  model  consists  of  parabolic  E  and  F  layers 

connected  by  a  linear  increase  in  electron  density  between  h _ E  and  a 

point  on  the  F^  parabolic  layer  where  the  plasma  frequency  is  equal  to 

I.  7  f  E.  This  model  is  described  in  detail  in  Table  II. 

o 

(9) 

Milsom'  1  showed  how  this  model  may  be  modified  slightly  using 
quasi  parabolic  E  and  F  layers  connected  by  a  "quasi  linear"  layer  of  the 
form 

N(r)  =  N  (r2  -  r2  )/ (r2  -  r2  )  (4- 1 ) 

m  b  m  b 

which  had  previously  been  shown  by  D.  B.  Muldrew^  and  D.  E.  Westover  ^  ^ 
to  allow  analytic  integration  of  the  ray  equations  in  terms  of  elementary 
functions. 

Note  that  the  Muldrew/ Westover /Milsom  quasi  linear  form  differs 
from  that  introduced  in  the  last  section.  Equation  (3.  12),  and  that  the  in¬ 
tegration  of  the  ray  equations  will  differ  as  a  consequence.  It  is  suggested 
that  the  functional  form  of  Equation  (3.  12)  is  also  a  suitable  candidate  for 
the  intermediate  layer,  with  the  advantage  that  solutions  to  the  ray  equations 
are  selected  from  a  single  subset  of  functions. 

In  order  to  compare  the  three  candidates,  the  Milsom  layer.  Equation  (4.  1) 
may  be  written  as 

N(r,  --  Ne  4  Ny  |rE/(rE  4  r  .!J  [(r/r,.)2  .  l]  (4.  2| 
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N(r)  =  N.  +  Nv  (1  -  r  /r)  (4.4) 

j  Y  Hj 

in  which  case  is  given  by, 

Ny  =  (N.  -  Ne)  •  r  /(r  -  r£)  •  (4.5) 

No  substantial  difference  obtains  between  intermediate  layers  (4.  2)  and 
(4.4);  we  do  present  the  specifications  based  upon  intermediate  layer  (4.4) 
in  Table  III;  the  graph  of  this  profile  is  shown  in  Figure  4.  1. 


TABLE  I. 

Definitions  of  Ionospheric  Terms 
r  -  radius  of  earth 

o 

f  E  -  E  layer  critical  frequency 

o 

h^E  -  E  layer  height 

Y  E  -  E  layer  semi-thickness 

m 

f  F,  -  F.  layer  critical  frequency 

h  F.  -  F.  layer  height 

ml  l 

v  F.  -  F.  layer  semi-thickness, 

7m  1  1  7  I 

f  F-,  -  F^  layer  critical  frequency 

h  F  -  F,  layer  height 

m2  2  7 

v  F,  -  F„  layer  semi-thickness 

7m  2  2 

N(r)  -  Electron  Density  Profile 

f  (r)  -  Plasma  Frequency  Profile 

kP  -  N(r)/f  (r)  =  1  .  2  4  E  4  (mhzf2  (cm) 

P 


J 


TABLE  II 


Specifications  for  Bradley-Dudeney  Profile 


N(r)  =  0 

r 

N(r)  =  Ne  -  Nx  (r/r£  -l)2 

VI 

x> 

u 

r 

N(r)  =  Ne  +  Ny  (r/rE  -1) 

rE  £ 

r 

£ 

N(r)  =  Nf  -  Nz  (r/rF  -l)2 

rj  * 

r 

S 

N(r)  =  0 

rt  £ 

r 

r  =  r  -  y  E 

b  E  ’  m 

r  =  r  +  h  E 

■£«  o  m 

rj  =  rF-VmF2U  -N./Nf)1/2 


r,„  =  r  -  y  F 
bF  F  7m  2 


Ne  =  k(foE)2 

nf  =  k(f/2)2 

N.  =  k  (I.  7  f  E)2 
J  0 

NX  =  NE<rE/>'mE)2 

NZ  =  NF<rF/>'mF>2 

NY  =  ,Nj  ‘  NE>  •  V'j  *  rE^ 
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TABLE  III 


Specifications  for  the  Modified  Bradley-Dudeney  Profile  (MBD) 


N(r) 

=  0 

N(r) 

=  ne 

-Nx(l 

-rE/r)2 

N(r) 

=  ne 

+  Ny(l 

-rE/r) 

N(r) 

=  nf 

-Nz(! 

-rF/r)2 

N(r) 

=  0 

£  rb 

£  rE 

£  r. 

J 


rj  =  rFrbF/{rbF+ymF2<1 
rF  '  r.  +  hmF2 

'«  =  rF  rbF/<rbF  *  ymF2> 

r„  =  r  -  y  F_ 

bF  F  m2 

=  k  (f  E)2 

NF  -  k<f„F2»2 

N.  =  k  (1.  7  f  E)2 
J  o 

NX  =  NE,rb/ymE)2 

NZ  '  NFfrbF/ymF2> 

Ny  =  (N.  -  Ne)  r./(r.  -  r£) 
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Modified  Bradley-Dudeney  Electron  Density  Profile 
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FIGURE  4-1 


5.  PIECEWISE  SMOOTH  QUASI  PARABOLIC  PROFILES 


A  Piecewise  Smooth  Quasi  Parabolic  Profile  (PSQP)  consists  of  alternating 
quasi-parabolic  segments  fitted  together  in  such  a  manner  that  the  electron 
density  and  its  first  derivative  is  continuous  at  each  fitting  point.  Consider 
two  such  segments  with  a  QP  segmented  fitted  below  an  inverted  QP  segment 


N(r)  =  Ne  -  Nb  (1  .  rE/r) 

r,  ^  r  £  r 
b  x 

N(r)  =  Nc  +  ND  (1  -  rc/r)2 

r  j;  r  ^  r 
x  *  *  y 

(5.  1) 


Continuity  of  the  electron  density  and  its  derivative  at  the  fitting  point  r 
demands  that 


NC  +  ND  ^  "  rc^rx)  =  NE  “  NB  ^  "  r£^rX* 


r  Nn  (r  -  r 
C  V  c  x 


=  rE  NB  (rx  "  rE  ) 


(5.2) 


or  equivalently. 


ND  =  NB  rE  (rx  '  V/rc  (rc  -  rx> 


NB  (rx  "  rE^  (rc  "  rEJ  =  (NE  "  Nc}  rc  rx 


(5.3) 


If  the  QP  segment  is  specified  ad  hoc  (N  ,  N  ,  r  are  given),  then  we 

E  13  E 

have  two  conditions  to  determine  four  parameters  (N_,  N„,  r  ,  and  r  ): 

C  u  c  x 

thus  two  parameters  may  be  determined  by  two  independant  arbitrary 
conditions.  For  example,  a  topside  (or  bottomside)  tail  maybe  constructed 
which  goes  smoothly  to  zero  (N£  =  0),  at  a  specified  point  (rc  is  given),  or 
fitted  at  a  specified  point  (r  is  given);  or  a  valley  may  be  constructed 
(Nc<  Na  is  given)  at  rc  or  fitted  at  r  .  The  algebra  to  solve  for  the 
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the  remaining  parameters  is  straight  forward. 


Between  two  QP  segments  a  valley  can  be  constructed  from  an  inverted 
QP  segment  fitted  smoothly  to  both  QP  segments.  In  this  case  we  have 
four  continuity  conditions  for  five  parameters  (N^.,  N^,  r^,  r^,  r  )  if  the 
top  QP  layer  is  given 


N(r)  =  Np  -  Ng(1  -  rp/r) 


r  £  r  <  r_ 

y  F 


To  determine  the  intermediate  layer,  one  may  arbitrarily  choose  the  valley 
depth  a  N  a  0  such  that  =  N^,  -  ^  N  and  the  remaining  parameters  are 
determined  from  the  four  continuity  conditions.  Again  the  algebra  is  straight 
forward.  Electron  density  profiles  with  *N  =0,  1/2  f  E,  and  f  E  are 

“  o  o 

shown  in  Figures  5.  1,  5.  2  and  5.  3. 

Another  interesting  criterior  to  impose  on  the  intermediate  layer  is  the 
scale  independent  condition,  expressed  by 

N  r  r  =  N  (r  r  )  (r  -  r  ) 

E  x  y  B  x  E'  '  y  E 

(5.5) 

N  r  r  =  N  (r  -  r  )  (r  -  r  ) 

x*  x  y  G  F  x  x1  y 

Note  that  each  of  these  two  conditions  are  independent  of  a  scale  change  of 
either  QP  layer;  the  solutions  for  the  fitting  points  and  r  depend  only  upon 
the  layer  heights  and  semi  thicknesses  of  the  two  QP  layers.  Thus,  in  an 
electron  density  model  where  only  f  E  and  f  F  vary,  the  fitting  points 
do  not  change.  The  algebra  is  straight  forward  to  demonstrate  that  these 
conditions  satisfy  the  continuity  conditions. 
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The  parameters  r^  and  r  are,  then,  given  by  the  solutions  to  the 

2 

quadratic  equation  Ar  -  B  r  +  C  =  D,  where 


rF 


E  B  '  E 


B  *  4  ng  (nb  -  ne»  -  4  NB  (ne  -  V 


[5.6) 


C  =  rE  rF  <rF  -  'E*  NB  NG 


An  electron  density  profile  satisfying  conditions  (5.5)  is  shown  in  Figure  5.4. 
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6. 


IONOSPHERIC  LAYER  MODELS  AND  ELECTRON 
DENSITY  PROFILES 


A  computer  program  has  been  written  which  performs  the  integration 
of  the  ray  equations  through  quasi-parobolic  ionospheric  layers.  The 
primary  source  of  models  for  the  individual  ionospheric  layers  used  in 
this  Fast  Ionospheric  Modeling  Program  (FIMP)  is  the  ITS  model  pro¬ 
viding  sets  of  numerical  coefficients  for  the  evaluation  of  foF£,  M^OOOJF^, 
11^2  and  f^E.  Each  ionospheric  characteristic  is  represented  for  a  given 
month  and  solar  activity  by  quasi  orthogonal-polynomial  expressions  in 
terms  of  geographic  latitude,  geographic  longitude,  and  coordinated  uni¬ 
versal  time.  Detailed  descriptions  of  this  model  are  found  in  the 
literature  1^,  15,  16,  17,  18^ 

This  model  of  ionospheric  characteristics  only  purports  to  provide 
monthly  median  values  of  the  individual  parameter;  on  any  given  day  devi¬ 
ations  up  to  20%  can  be  expected.  The  FIMP  program  expects  input  data 
taken  from  a  recent  vertical  ionogram  from  which  the  model  character¬ 
istics  may  be  updated  to  more  nearly  reflect  the  actual  ionospheric  con¬ 
ditions  prevailing  at  the  time  for  which  the  simulation  study  is  being 
constructed. 

The  updating  process  employed  is  quite  simple:  if  Pq  is  the  observed 
value  of  a  given  ionospheric  parameter  extracted  from  the  vertical  ionogram, 
and  P^  is  the  model  value  of  that  parameter  at  the  time  and  place  of  the 

vertical  ionogram,  then  a  correction  factor  k  =  P  / P  is  evaluated. 

6  '  pom 

Subsequently,  when  that  parameter  P^  is  extracted  from  the  model  at 
some  other  (nearby)  time  and  place,  the  parameter  is  updated  to  the 
value  given  by  P  =  k  P  . 

°  7  n  n  m 


At  the  beginning  of  a  run  of  the  FIMP  program,  information  from  the 
vertical  and  oblique  ionograms  are  input,  and  a  table  of  required  ionospheric 
characteristics  is  constructed  along  the  great  circle  defined  by  the  trans¬ 
mitter  and  its  boresight  azimuth. 

Information  from  the  oblique  ionogram  may  also  be  used  to  correct 

the  ionospheric  characteristics  model,  in  the  sense  of  constructing  a  down 

range  model  of  f  F_  constrained  to  reproduce  the  leading  edge  of  the 
o  ^ 

oblique  ionogram. 

This  technique  utilizes  the  initial  climatological  variation  of  all  para¬ 
meters  except  for  f  F_.  Information  from  the  oblique  ionogram  is  input 
to  the  FIMP  program  as  a  list  of  (one  way)  group  paths  corresponding  to 
a  set  of  equally  spaced  frequencies.  Initially  f qF ^  is  established  by  the 
updated  model  as  well  as  the  other  parameters  of  the  synthetic  ionospheric 
profile.  The  first  frequency  in  the  set  is  selected. 

The  minimum  group  path  corresponding  to  that  frequency  fQ  in  that 
synthetic  profile  characterized  by  fQF ^  is  calculated  and  compared 

with  the  corresponding  group  path  of  the  leading  edge.  In  general  these 
two  values  would  not  agree;  assume  that  the  simulated  minimum  group 
path  is  larger  than  the  leading  edge  group  path. 

The  simulated  minimum  group  path,  in  this  case,  would  be  appropri¬ 
ate  on  the  actual  oblique  ionogram  for  a  larger  frequency  f  than  the 
operating  frequency  fQ.  This  frequency  f  may  be  determined  from  inter¬ 
polation  within  the  array  of  leading  edge  group  paths. 

The  discrepancy  between  simulated  and  observed  group  paths  is  now 

ascribed  to  an  error  in  the  value  of  f  in  the  synthetic  profile;  in  this 

c 
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example  f  is  too  small  and  should  be  increased.  A  new  value  of  critical 
frequency  is  then  defined  by 


and  the  process  is  iterated  until  the  discrepancy  between  simulated  and 
observed  group  paths  is  small  (less  than  0.  5  km). 

The  limit  of  this  iterative  process  results  in  a  value  of  the  critical 
frequency  f  which  is  then  presumed  to  be  the  appropriate  value  of  f^F^ 
at  the  reflection  point  of  the  simulated  minimum  group  path.  As  the  fre¬ 
quencies  in  the  defined  set  are  so  processed,  a  function  f  F_  versus  range 
from  transmitter  is  developed  which  reproduces  the  leading  edge  of  the 
oblique  ionogram  within  the  established  tolerance. 

However,  this  function,  despite  its  ability  to  reproduce  the  leading 
edge  so  well,  is  not  considered  reliable  in  detail,  but  is  considered  reli¬ 
able  to  provide  the  down  range  trends  of  f^F^.  ^he  results  of  the  previous 
analysis  are  fitted  to  a  straight  line  constrained  to  fit  the  updated  model 
at  the  transmitter  location.  Furthermore,  only  elements  of  the  function 
at  ranges  exceeding  1000  kilometers  are  employed  in  the  fitting  procedure. 
This  linear  model  for  f^F^  is  then  employed  in  a  simulation  of  the  leading 
edge,  discrepancies  from  the  observed  ionogram  are  taken  as  reliability 
estimates  of  the  FIMP  system.  The  linear  model  is  further  employed  in 
further  range  registration  calculations  desired  by  the  user  for  the  current 
program  run. 

The  electron  density  profile  chosen  for  the  FIMP  ray  calculation  pro¬ 
gram  is  a  piecewise  smooth  quasiparabolic  model  based  upon  f^E,  ^F^, 

h  F,,  and  y  F  being  modeled  from  the  global  ionospheric  characteristic 
m  £  me. 
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model.  The  layers  consist  of  (1)  an  inverted  QP  ”D"  layer  attached  to  the 
E  layer  at  the  100  km  level  falling  smoothly  to  zero  below  that  level, 

(2)  a  QP  E-layer  at  a  layer  height  100  km  and  semithickness  20  km, 

(3)  an  inverted  QP  layer  attached  between  the  E  and  F  layer  according 
to  the  scale  independent  criterior  established  in  Section  5,  and  (4)  a  QP 
F  layer.  In  addition,  the  ray  on  the  way  down  may  encounter  a  different 
E-layer  from  that  on  the  way  up. 

Since  each  ray  in  the  frequency-elevation  angle  family  probes  in  gen¬ 
eral  a  different  region  of  the  ionosphere  each  such  ray  requires  a  different 
profile.  The  profile,  then,  is  synthesized  from  fQE  for  the  upward  track 
taken  at  the  point  of  joining  the  intermediate  layer,  f^E  for  the  downward 
track  similarly  chosen,  and  the  F  layer  parameters  chosen  at  the  reflec¬ 
tion  point.  Because  the  solution  for  these  points  are  known  a  posteriori, 
the  final  solution  is  approached  in  an  iterative  fashion.  This  iteration 
does  converge  rapidly,  especially  if  care  is  exercised  to  select  a  frequency- 
elevation  member  of  the  ray  family  close  to  the  most  recently  completed 
ray  calculation. 


7.  APPLICATIONS  TO  TRANSMISSION  CIRCUITS 


One  of  the  most  common  applications  of  ionospheric  ray  tracing  pro¬ 
grams  is  the  determination  of  the  circuit  parameters  (transmission  fre¬ 
quency,  f  and  elevation  angle,  j3)  required  for  communic'ition  between  given 
transmitter  and  receiver  locations.  For  an  arbitrarily  specified  ionosphere 
this  usually  involves  the  use  of  a  homing  algorithm  with  a  three  dimensional 
ray  tracing  program.  When  the  electron  density  gradients  are  large  even 
this  technique  does  not  yield  unique  results  and  the  computer  time  required 
to  find  any  possible  circuits  may  be  quite  large.  If  an  analytic  expression 
for  the  range  in  terms  of  f  and  8  can  be  used  the  problem  may  be  greatly 
simplified.  The  quasi-parabolic  model  has  been  used  to  obtain  solutions  to 

the  homing  problem  for  several  ionospheric  models.  The  results  we  pre¬ 
sented  include  those  with  several  layers  but  neglect  the  effects  of  the 
horizontal  gradients. 


The  simplest  case  is  that  of  a  single  layer  (SLOP  night  time  iono¬ 
sphere)  when  the  ionosphere  is  entirely  specified  by  fQF  ,  h^F^  ant*  ^m^ 2’ 
We  characterize  the  transmission  curves  by  the  ratio  of  the  transmission 
frequency  f,  to  :  -  ^^0^2* 

The  ground  is  given  by  the  following  expression. 


2a  |  (Y-P)  -  F  a.cos_llogo  /Ur 
'  2V  C  UV2 


where 


r  r»  r 

b  m 1 


„2  2  2 

-Fa  cos  p 


m 
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A  =  F  -  1  +J  b  | 

\m) 

B  -  -2  r  /rb  \2 

”(W 

U  =  B2  -  4  AC 


W  =  2\ZcT  F  r,  sin  Y+  2C  +  B  r, 
b  b 


cos  Y  = -  cos  3  ,  Y  is  the  elevation 

rb 

angle  at  the  base  of  the  ionosphere 


r  =  a  +  h  F_ 
m  m2 


r,  =  r  -  y  F-  and  a  is  the  radius  of  the  earth, 
b  m  m  2 


The  group  path,  P'  is  given  by 

P'  =  2  |l  -  F2/a|  rb  siny  -2a  sing 
where  V  =2  Ar^  +  B  +  2  r^  F\T  A  sin  Y  . 


Figures  7-1  and  7-2  show  the  ground  range  and  group  path  for  a  sample 

ionosphere  with  h  F_  =  325  km.  and  y  F_  =  175  km.  The  elevation 
r  m2  m  2 

angles  are  sampled  at  1°  intervals  for  F  values  between  0.2  and  3.0  at 
intervals  of  0.2.  The  determination  of  the  circuit  parameters  for  a 
specified  range  Rq  consists  in  finding  the  solutions  of  the  equation 

Ro  =  R(F,  8  ) 

We  first  find  the  solution  for  the  smallest  possible  elevation  angle.  There 
are  two  cases: 

i)  Specular  reflection  from  the  bottom  of  the  ionosphere, 

F  =  0+  P  .  .  tan'1  MV4  3/rbl 
m,n  - ain/R  72 a\ - 


:  /2a\ 
°  ) 
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ii)  F  4  0  B  .  =  0°  . 


mm 


Other  solutions  are  then  searched  for  at  elevation  angles  greater  than  3 


min 


for  a  given  angular  step  size.  For  each  new  elevation  angle  an  iterative 
procedure  is  used  starting  with  the  value  of  F  found  as  a  solution  for  the 
previous  angle.  A  correction  to  this  F  is  found  using: 


(R(F  B)  -  R  ) 

AF.  = - »  - 5_ 

1  3  R 


■f^(Fi.B) 


F.  +  l  =  F.  +  AF. 

i  ii 


where  3R  is  calculated  analytically  using 

ST 


where 


3R 

3F 


cosB  jl 


Zy/ZT 
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2C  3F  " 


1 
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3U 

dF 
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3  A  ,  _  and 

ST‘  2F 


3  W  =  2y/c~ 
3F 


r^  sinY+|F  r^  siny+  2] 


VC 


3  C 
3F 


We  search  for  solutions  for  a  given  range  in  l/4°  steps  in  elevation  angle 
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between  8  .  and  8  ,  the  angle  at  which  the  ray  penetrates  the  .maximum 

min  max  °  7  r 

height  of  the  layer.  This  occurs  where  U  5  0  which  yields 


For  this  sampling  internal  the  procedure  converges  rapidly  using  the  start¬ 
ing  frequency  of  the  previous  angle.  For  a  given  range  the  solutions  include 
the  low  rays,  the  high  rays  and  the  MUF.  Figure  7-3  shows  all  the  solutions 
obtained  for  a  3000  km  path  using  the  sample  set  of  ionospheric  parameters 
whose  transmission  curves  where  shown  in  Figure  7-1.  We  have  considered 
propagation  by  up  to  four  hop  ionospheric  modes.  The  first  hop  starts  from 
initial  condition  ii),  while  the  solutions  for  the  other  three  hops  were  ob¬ 
tained  beginning  with  initial  condition  i).  For  each  pair  of  circuit  para¬ 
meters  on  the  curves  we  calculate  the  group  path.  The  resulting  bistatic 
ionogram  is  plotted  as  a  function  of  the  normalized  frequency,  F  in 
Figure  7-4.  In  a  particular  application  one  would  only  consider  that  part 
of  the  ionogram  that  would  be  possible  for  the  given  transmitter  antenna 
pattern.  While  we  do  not  explicitly  account  for  the  effects  of  horizontal 
gradients  in  the  ionogram  it  is  possible  to  obtain  some  indication  of  the 
size  of  these  effects  by  considering  other  stratified  ionospheric  layers 
with  parameters  that  are  slightly  different  from  the  original  problem. 

For  the  test  ionosphere  this  might  include  a  range  of  300-350  km  for 

h  F  and  150-200  km  for  y  F,. 
m  ^  m  <- 

When  we  use  a  multi-layer  quasi -parabolic  model  (MLQP)  we  have 
in  addition  to  the  layer  parameters,  three  other  parameters  for  the 
underlying  ionization  of  the  F  layer.  We  have  taken  constant  values  of 
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FREQUENCY  <f> 


FIGURE  7-3 

a 


ELEVATION  ANGLE  38 


GRQ 


110  km  and  20  km  for  the  height,  h  E  and  semithickness,  y  E  of  the  E  layer 

0  m  m 

and  examine  the  effect  of  various  values  of  f  E  upon  the  transmission  curves. 

o 

In  Figure  7-5  we  show  the  vertical  electron  density  distribution  obtained 
for  various  combinations  of  the  layer  parameters.  We  have  fixed 


F? = 
m  2 

325  km  for  all 

six  profiles. 

The  other  parameters  used 

f  F  _ 
o  2 

y  F? 

m  2 

f  E 
o 

(a) 

10.  0MHz 

175.  km 

3.  5MHz 

(b) 

10.  0 

125. 

3.  5 

(c) 

10.  0 

100. 

3.  5 

•  r - —  —  • 

(d) 

10.  0 

175. 

2.5 

(e) 

7.  0 

175. 

3.  5 

(f) 

7.  0 

125. 

3.  5 

We  have  calculated  the  range  and  group  path  for  frequency  intervals  of 
2MHz  and  angular  intervals  of  1°  for  the  parameter  sets  (a)-(e). 

Figures  7-6  and  7-7  are  plots  of  R  and  P1  for  the  parameters  in  (a) 
and  Figures  7-8  through  7-10  are  the  range  curves  for  parametrs  (b), 

(d)  and  (e).  The  maximum  frequency  that  can  propagate  by  the  E  layer 
is  given  by: 

For  case  (d)  f  E  =  2.  5MHz  and  f  E  =  1  3.  4MHz  and  for  all  the  other  cases 
o  max 

f  E  =  3.  5MHz  and  f  E=18.  8MHz.  The  absence  of  E  modes  above  these 
o  max 

frequencies  is  clearly  visible  in  Figures  7-6  through  7-10.  The  E  and  F 
propagation  modes  are  clearly  separated  in  all  of  the  plots.  With  a  1° 


f  E  = 
max 
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sampling  interval  in  |3we  find  only  one  or  two  high  angle  E  rays  for  each 
set  of  parameter  for  any  fixed  frequency  of  transmission.  A  comparison 
of  Figure  7-1  with  7-6  reveals  that  the  E  layer  has  very  little  effect  on  the 
range  curves  for  most  of  the  high  rays  as  we  get  away  from  the  skip  dis¬ 
tance.  For  a  fixed  frequency  we  also  see  from  Figure  7-6  and  7-7  that 
the  skip  distance  and  minimum  group  path  do  not  occur  at  the  same  eleva¬ 
tion  angle. 

From  Figure  7-3  we  see  that  for  the  SLQP  model  if  there  is  a  solution 
to  the  equation  RQ=R(f,  (3)  that  it  is  a  unique  solution.  In  the  case  of  the 
MLQP  model  the  crossing  of  the  various  frequency  curves  implies  that 
it  is  possible  to  have  dual  frequency  solutions  for  the  same  elevation  angle 
for  F  mode  propagation.  This  complication  has  lead  us  to  use  a  different 
method  for  finding  all  the  possible  solutions  for  a  given  range  problem.  In 
this  case  we  attempt  to  find  the  elevation  angle  that  completes  the  path  for 
a  fixed  frequency  rather  than  searching  for  the  frequency  for  a  fixed  eleva¬ 
tion  angle.  For  each  frequency  we  use  the  electron  density  profile  to 
compute  the  limits  on  the  initial  elevation  angles  for  reflection  from  the 
different  layers.  This  results  in  a  set  of  angles: 


minimum  elevation  angle  for  E  layer  reflection. 

elevation  angle  for  E  layer  skip  distance. 

maximum  elevation  angle  for  E  layer  reflection 
or  minimum  elevation  angle  for  F  layer  reflection 

elevation  angle  for  F  layer  skip  distance. 

maximum  elevation  angle  for  F  layer  reflection 
(penetration  angle). 


Between  any  successive  pair  of  these  angles  the  range  curve  is  a  monotonic 


function.  The  range  at  the  end  points  determines  whether  a  solution  exists. 
A  search  procedure  is  then  initiated  to  find  the  elevation  angle  which  results 
in  a  range  within  a  small  error  of  the  specified  circuit  range. 


A  single  method  is  used  to  find  the  appropriate  solution  even  though 
has  a  wide  range  of  variation.  The  calculation  is  per- 


the  function 


formed  with  the  variable  cos  3  rather  than  3  itself.  The  ground  range  is 

2 

calculated  for  a  cos  3j  midway  between  the  values  at  the  end  points.  De- 


2 

pending  upon  the  sign  of  (R  -R(f  ,  cos  {3^)  )  one  of  the  end  points  is  replaced 
2 

with  cos  3^  and  a  new  midpoint  is  taken.  The  procedure  is  terminated 
when 


|RQ-R(f,  cos2^^)  (  <0.5  km. 

Figure  7-11  is  a  plot  of  F  mode  solutions  for  Rq=  3000  km  for  ionospheric 
model  (a).  Figure  7-12  is  a  plot  of  the  group  path  as  a  function  of  frequency 
for  these  solutions.  These  curves  corresponds  to  the  same  ionosphere 
parameters  that  were  used  in  the  SLQP  model  with  the  addition  of  an  under- 
lying  E  layer.  Comparing  Figure  7-11  with  Figure  7-3  we  see  that  most 
of  circuit  parameters  for  rays  above  the  MUF  for  each  hop  coincide  with 
each  other  for  the  two  ionospheric  models.  For  the  MLQP  model  we  only 
searched  for  solutions  for  f  2  {  F_.  The  effect  of  changing  individual 
model  parameters  upon  the  propagation  parameters  and  the  resulting 
bistatic  ionogram  for  the  F  modes  is  shown  on  Figures  7-13  through  7-16. 

In  Figures  7-13  and  7-14  only  the  semithickness  of  the  F^  layer  has  been 
changed  (case  (b))  while  in  Figures  7-15  and  7 - 1 6  (case  (d))  the  value  of 
foE  has  been  changed  from  3.  5MHz  to  2.  5MHz.  Both  changes  have  very 
little  effect  on  the  propagation  above  the  F  layer  MUF's.  Reducing  the 
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semithickness  results  in  raising  the  elevation  angle  necessary  to  reach  the 
receiver  site  for  a  given  frequency  for  the  low  angle  rayp.  This  causes  an 
increase  of  the  group  path  of  ~50  km  for  the  one  hop  mode  and  a  correspond¬ 
ing  increase  for  each  of  the  multihop  modes.  Reducing  the  value  of  f  £ 
allows  for  F  mode  propagation  at  smaller  elevation  angles  to  the  receiver. 
For  these  elevation  angles  the  group  path  is  reduced  by  50-100  km,  inde¬ 
pendent  of  the  number  of  ground  hops. 

An  examination  of  the  range  curves  for  models  (a)  -  (d)  shows  that  prop¬ 
agation  over  distances  of  1000-3000  km  are  characterized  by  F  layer  MUF's. 
By  reducing  the  value  of  f^F^  we  can  observe  the  occurence  of  E  layer 
MUF's.  We  have  used  case  (e)  with  a  ground  range  of  2200  km  to  display 
these  modes.  Figure  7-17  shows  the  circuit  parameters  required  for  this 
range  for  both  the  one  and  two  hop  E  and  F  modes.  Figure  7-18  is  the  result¬ 
ing  group  path  for  these  parameters.  The  four  individual  MUF's  are  shown 
on  Figure  7-17.  They  are: 

One  Hop:  F  MUF  =  17.0MHz. 

E  MUF  =  18.  3 
Two  Hop:  F  MUF  =  10.  8 

E  MUF  =  13.  0 

The  method  of  searching  for  solutions  to  the  cirauit  problem  is  suffi¬ 
ciently  sensitive  so  that  all  mathematically  possible  solutions  are  obtained 
including  those  which  yield  very  little  energy  at  the  receiver,  i.  e.  when 

dR  is  very  large.  Some  of  the  solutions  shown  on  Figure  7-17  are  not 
b cos 2 8 

apparent  when  the  range  curve  of  Figure  7-10  is  examined.  An  example 
of  this  is  the  lowest  frequency  low  angle  1  F  mode  which  almost  coincides 
with  a  high  angle  2  E  mode  between  14°  and  15°  just  under  12MHz. 
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The  MLiQP  model  has  been  used  to  simulate  a  discrete  oblique  ionogram 
for  very  long  ranges  in  conjunction  with  a  study  of  high  frequency  ducted 
ionospheric  modes.  Ray  tracing  programs  with  electron  density  models 
that  include  horizontal  gradients  have  been  used  to  determine  ionospheric 
conditions  that  are  favorable  for  ducted  propagation.  The  MLQP  model 
has  been  used  to  study  the  effect  of  multihop  ground  modes  on  the  discrete 
oblique  ionogram.  Figure  7-19  shows  the  circuit  parameters  for  F  mode 
multihop  propagation  over  a  ground  range  of  14,000  km.  Figure  7-20  is  a 
plot  of  the  group  path  resulting  from  these  parameters.  The  calculations 
were  performed  using  the  layer  parameters  of  model  (a). 
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FIGURE  7-19 
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8.  SIMULATION  OF  VERTICAL  AND  OBLIQUE 
IONOGRAMS _ _ 

Penultimately  the  objective  of  creating  ionospheric  electron  density 
models  and  rapid  ray  calcuation  procedures  is  to  provide  accurate  simu¬ 
lation  of  RF  propagation  characteristics  in  the  ionosphere  prevailing  at 
the  time  and  place  of  the  simulation.  If  this  accuracy  objective  is  met, 
the  model  may  then  be  employed  as  a  reliable  tool  to  estimate  system 
performace  characteristics,  and  target  range.  This  section  describes 
how  the  electron  density  model  may  be  finely  adjusted  based  upon  the 
gross  and  detailed  features  of  vertical  and  oblique  ionograms. 

The  gross  features  of  the  vertical  ionogram  are  used  to  correct  the 
ionospheric  parameters  fQE,  fQF2,  ^'F2'  M(^000)F2  *rom  which  a  profile 
is  constructed.  The  vertical  ionogram,  virtual  height  (group  path)  versus 
frequency,  is  then  generated. 

Five  simulations  of  the  vertical  ionogram  are  shown  in  Figure  8-1. 

In  all  of  these  simulations  the  primary  ionospheric  parameters  are  as 
follows: 


f  E  = 
0 

3.  8Mhz 

f  F  = 
0  Z 

13.  6Mhz 

h  E  = 
m 

1  1 0.  0km 

h  F  = 
m  Z 

333,  0km 

y  e  = 

7m 

20.  0km 

y  F0  = 
m  2 

1  06.  4km 

The  differences  lie  in  the  intermediate  layer  model  between  the 
E  and  F  layer.  No  ionization  produces  ionogram  A.  Ionogram  B  results 
from  a  constant  electron  density  equal  to  that  at  the  E  layer  height  up  to 
the  intersection  with  the  F  layer.  Ionogram  C  is  produced  from  an  in¬ 
verted  QP  layer  joining  the  E  layer  smoothly  at  h^E  1°  the  intersection 
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on  the  F  layer  where  the  plasma  frequency  equals  1.  7  f  E,  while  C  results 

o 

from  the  modified  Bradley-Dudeney  profile.  The  piecewise  smooth  QP 
layer  with  the  scale  independant  criterior  results  in  ionogram  D.  The  fit¬ 
ting  points  and  plasma  frequencies  for  these  models  are  listed  in  Table  8.  1. 

It  is  clear  from  inspection  of  Figure  8-1  that  the  choice  of  electron 
density  model  between  the  E  and  F  layer  has  a  considerable  effect  upon 
the  detailed  structure  of  the  vertical  ionogram. 

The  piecewise  smooth  QP  profile  using  the  scale  independant  criterior 
(ionogram  D)  tends  to  flatten  out  the  belly  from  the  QP  layer  model  (iono¬ 
gram  A),  while  a  shelf  tends  to  strongly  lift  the  curve  just  above  f^E 
(ionogram  B).  Electron  density  Model  C  is  close  to  the  Bradley-Dudeney 
profile  except  that  the  intermediate  layer  is  tangent  to  the  E  layer  at 
h^E;  ionogram  C  tends  to  rise  more  or  less  linearly  just  above  fQE, 
characteristic  of  the  inverted  QP  behavior  to  a  cusp  singularity  at  f.. 

Above  F  the  ionogram  displays  the  concavity  characteristic  of  the  QP 
layer.  Note,  however,  that  the  cusp  at  f.  is  much  less  severe  than  that 
produced  by  the  Bradley-Dudeney  profile. 

The  scale  independent  criterion  for  the  construction  of  the  intermedi¬ 
ate  layer  generates  a  fitting  point  high  up  on  the  F  layer,  while  the  iono¬ 
gram  to  be  modeled  displays  a  pronounced  belly  similar  to  ionogram  A. 
Adjustment  of  the  F  layer  height  and  semi  thickness  does  not  overcome 
this  difficulty,  which  reappears  in  the  context  of  simulating  the  oblique 
ionogram,  discussed  below. 

Clearly,  the  question  of  the  intermediate  ionization  model  requires 
more  study  using  a  larger  sample  of  ionograms,  and  the  systematics  of 
selecting  the  best  intermediate  layer  model  developed. 


TABLE  8-1 


Parameters  for  Intermediate  Layer  Models 


Model 

Layer 

hI 

hL 

hJ 

fI 

f 

L 

£J 

A 

130 

130 

227 

0 

0 

0 

B 

1 10 

110 

231 

3.  80 

3.  80 

3.  80 

C 

1 10 

110 

239 

3.  80 

3.  80 

6.  46 

C' 

110 

- 

239 

3.  80 

- 

6.46 

D 

112 

127 

282 

3.  78 

3.  60 

11.  98 

The  main  feature  of  the  oblique  backscatter  ionogram  to  simulate  is 
the  leading  edge.  The  basic  problem  is,  for  a  given  ionosphere  and  at  a 
given  frequency,  find  that  ray  (selected  by  elevation  angle)  that  corresponds 
to  the  minimum  group  path  length  (MGP),  or,  equivalently,  locate  a  zero 
in  the  derivative  of  the  group  path  function  of  take  off  elevation  (or  its 
cosine).  We  shall  discuss  the  MGP  solution  for  F  layer  reflections;  there 
would  be  no  substantial  difference,  however,  should  one  desire  the  E 
layer  MGP  solutions. 

A  minimum  in  the  group  path  length  occurs  because  of  the  maximum 
in  the  electron  density  at  the  F  layer  height.  As  one  increases  the  ele¬ 
vation  angle  from  low  angles  the  group  path  decreases,  and  to  a  first  ap¬ 
proximation  the  situation  is  similar  to  single  reflection  from  a  fixed  plane 
mirror.  This  simple  model  is  inadequate  in  two  respects.  Firstly,  the 
actual  penetration  into  the  ionosphere  increases  with  elevation,  and 
secondly,  the  time  delays  increase.  Thus,  we  need  to  float  the  mirror 
upwards  as  we  increase  the  elevation  to  keep  pace  with  the  virtual  height, 
tending  to  increase  the  group  path  length.  As  the  virtual  frequency 
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approaches  the  critical  frequency  of  the  layer,  the  group  path  length  in¬ 
creases  without  bounds;  this  increase,  then,  will  eventually  dominate  the 
previous  decrease  and  thus  a  minimum  in  the  group  path  length  will  occur. 
Rays  at  elevations  less  than  that  corresponding  to  the  minimum  are  the 
"low  angle"  rays,  for  which  a  geometric  path  loss  model  applies  quite 
well.  The  group  path-elevation  curve  increases  quite  sharply  above  that 
angle,  resulting  in  rapid  decrease  in  energy  density  and  additional  losses 
must  be  allocated.  Note  that,  at  the  minimum,  focusing  occurs  resulting 
in  signal  enhancement. 

During  traversal  of  a  ray  through  an  inverted  QP  layer,  the  group 
path  decreases  with  increasing  elevation  angle,  it  is  only  after  entering  a 
QP  layer  when  the  group  path  increases.  This  establishes  a  search  window 
for  the  minimum  group  path  solution  between  penetration  into  the  layer 
and  penetration  through  the  layer. 

Penetration  into  the  F,  layer  occurs  when  the  function  Z(r  )  given  by 
quation  3.4  vanishes,  where  r.  is  the  point  at  which  the  intermediate 
inverted  QP  layer  meets  the  QP  F~  layer,  or  at  the  angle  B.  given  by 

Lt  OjJ 

cos2g  .  =  (r.  Ir  )2  (1  .  kN(r  .)/f2) 
oj  Jo  J 

Penetration  through  the  F^  layer  occurs  when  the  discriminant  of  the 
function  Z(r  )  given  in  equation  3.4  vanishes,  or  at  the  angle  P  given  by 

c°S\f  =  XG  "  -  V'"  -  XC  +  XG' 

where  the  other  terms  are  as  defined  following  equation  3.  4. 

The  convenient  independant  variable  to  employ  in  lieu  of  the  takeoff 

elevation  anele  B  is  the  square  of  its  cosine,  x  -  cos^P  .  The  solution 
^  o  o 


for  the  minimum  group  path  at  x  will  be  found  in  the  domain  given 


by 


X.  2  X  2  X 

J  rn  F 


This  interval  is  divided,  by  convenience,  into  eight  equal  interval  seg¬ 
ments  of  width  A  =  (Xj_.  -  Xj)/8.  To  start  the  search  process,  three  rays 

calculated  yielding  group  paths  P  =  P( x  )  where  x  =  x.  +  kA,  k  =  1,  2,  3. 

k  K  k  J 

Numerical  estimates  of  the  derivatives  of  P(x)  are  calculated  as 


are 


Pk  = 

P’(XA>  =  (P2-P1)/(X2‘ 

V 

PB  = 

p,(xB)  =  (P3  -  P2)/(x3  - 

X2) 

where 

X  A  = 

(X1  +  x2)/2 

and 

X  B  = 

(  x 2  +  x3)/2 

The  first  and  second  derivatives  of  the  function  P(x)  are  estimated  by 

=  P'(V  =  (PX  +  PJ,)/2 

p„  .  p,,Xo)  =  ,pk  -  Pis)/(*A  -  «B) 

where  x  =  (x  +  x  )/2. 

o  O 

In  the  neighborhood  of  x  the  group  path  function  P(x)  may  be  ap¬ 
proximated  as 

P(X)  =  P(X  )  +  P'(  x  )(x_x  )  4  1/2  P"(X  )  ( x  _  x  )2 
o  o  o  o  o 

and  its  derivative  by 

P'(x)  =  P'(*D)  +  P"(xo)(x  -  Xq) 

At  the  minimum  of  P(x)  its  first  derivative  must  vanish,  thus  a  cor¬ 
rection  to  x  is  generated  as 

6x  r  -  p«  IP" 


6b 


The  points  x-,  and  x are  demoted  to  the  roles  of  points  x  and  x  respec - 

“  j  1  z 

tively,  while  a  new  value  of  is  given  by 

x„  =  x  +  6  x. 

3  o 

A  new  ray  is  calculated  yielding 

The  iteration  continues  until  is  within  0.  5  kilometers  of  the  most 
recently  calculated  minimum  value,  except  that  no  more  than  eight  rays 
are  calculated.  Escape  from  the  iterative  process  usually  occurs  after 
the  minimum  of  four  rays;  a  single  correction  is  generally  sufficient. 
Rarely  is  it  necessary  to  calculate  six  or  more  rays. 

Having  located  the  minimum  group  path  solution  at  the  given  fre¬ 
quency  for  the  specified  model  ionosphere,  the  model  ionosphere  may  be 
updated  to  take  into  account  horizontal  variations  of  the  ionospheric 
characteristics,  both  by  creating  a  new  synthetic  ionosphere  and  by 
calculating  layer  by  layer  changes  in  the  elevation  angle  due  to  the  hori¬ 
zontal  gradients.  With  this  new  ionosphere,  the  minimum  group  path  al¬ 
gorithm  may  be  reentered  to  improve  the  soltuion.  This  process  is 
reiterated  until  the  desired  convergence  criterior  is  attained. 

This  iterative  ionospheric  parameter  updating /direction  cosine  cor¬ 
rection  process  could  be  accomplished  within  the  iterative  search  for 
minimum  group  path  ray  by  ray  or  even  layer  by  layer  within  a  ray  cal¬ 
culation.  This  may  speed  up  the  iteration  process;  on  the  other  hand, 
such  a  double  correction  may  lead  to  instabilities  or  oscillations,  and 
may  actually  slow  up  the  process.  Furthermore,  the  parameter  correction 
scheme  to  fit  a  given  oblique  ionogram  depends  upon  the  result  of  the 
simulation  of  the  leading  edge  point,  and  the  double  iteration  is  required 


The  technique  described  in  Section  6  for  the  development  of  a  f 
versus  range  function  was  applied  to  the  simulation  of  the  leading  edge  of 
an  actual  oblique  ionogram.  The  function  so  developed  is  shown  in  Figure  8.  2 
together  with  the  ITS  model  for  foF-,  (corrected  to  the  vertical  ionogram 
parameters).  Also  the  linear  fit  f^F^  is  displayed. 

The  forced  fit  model  exhibits  oscillations.  This  behavior  may  in  fact 
be  real,  due  perhaps  to  an  ionospheric  disturbance  with  a  500  kilometer 
wavelength.  On  the  otherhand,  the  oscillations  may  be  due  to  compu¬ 
tational  prejudices  such  as  arbitrarily  selecting  f  F_  as  the  only  para- 

o  ^ 

meter  allowed  to  deviate  from  the  ITS  model.  Clearly,  h  F_  and  y  F_ 

7'  m  2  m  2 

have  significant  influence  in  determining  the  shape  of  the  leading  edge. 

The  oblique  backscatter  ionogram  was  simulated  using  three  models 
for  (1)  the  ITS  climatological  model,  (2)  the  linear  fit  model  based 

upon  the  forced  fit  process,  and  (3)  the  overhead  value  of  ^F^  at  the 
transmit  site.  The  results  of  these  simulations  are  shown  in  Figure  8.  3, 
and  residuals  are  given  in  Table  8.  2.  Clearly,  the  linear  fit  model  is 
significantly  better  than  the  ITS  model.  The  transmit  site  value  pro¬ 
duces  results  nearly  as  good  as  the  linear  fit  model;  however,  this  may 
only  be  fortuitous. 
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TABLE  8-2 


Residuals  of  Oblique  Ionogram  Simulation 

With  Various  f  F~  Models 
o  2 


Frequency  Residuals  From 


MH 

z 

Overhead  Value 

Linear  Fit  Model 

ITS  Model 

14 

+  7.8  km 

+  3.0  km 

+  15.8  km 

16 

+47.  6 

+42.  3 

+  65.3 

18 

+42.  0 

+36.2 

+  71.2 

20 

+21.8 

+  15.  5 

+  65.2 

22 

-  1.7 

-  8.6 

+  58.0 

24 

+  12.  6 

+  4.4 

+  87.8 

26 

+36.6 

+26.6 

+  131.  3 

28 

+72.4 

+60.  2 

+  192.  1 

30 

+72.  8 

+58.  0 

+226. 0 

32 

+41.8 

+23.  8 

+241. 7 

34 

+35.  2 

+  13.  2 

+289. 1 

36 

-38.  1 

-68.  2 

+295.2 

38 

-12.  7 

-55.  9 

+434. 3 

GROUP  PATH  (km) 
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